La solución más conocida para este problema es $x= e^{-\lambda}x_0$. Cuya gráfica es

import matplotlib.pyplot as plt
import numpy as np
plt.style.use('https://github.com/dhaitz/matplotlib-stylesheets/raw/master/pitayasmoothie-dark.mplstyle')
lamda = 5
x0 = 10
t = np.arange(0,10,0.01)
solrealy = np.exp(-lamda*t)*x0
plt.plot(t,solrealy)

[<matplotlib.lines.Line2D at 0x1876e5a8d50>]

def EF_caso1(lamda,h,N,x0):
    return ((1-lamda*h)**N)*x0
def EB_caso1(lamda,h,N,x0):
    return (x0*(1/(1+lamda*h))**N)
def CN_caso1(lamda,h,N,x0):
    return (x0*((2-lamda*h)/(2+lamda*h))**N)
hs = [0.1,0.2,0.03,0.4,0.5,0.75,1]
ts = [np.arange(0,10,i) for i in hs]

lamda = 5
x0 = 10

mat_parte1 = np.ones([3,len(hs)]) # 3 metodos* n hs * n tiempo
mat_parte1 = [list(mat_parte1[i]) for i in range(3)]
metodos_num_caso1 = [EF_caso1,EB_caso1,CN_caso1]
metodos = ['EF','EB','CN']
num = 0
plt.figure()
fig,axs = plt.subplots(1,3,layout = 'constrained')
for metodo,ax in zip(metodos_num_caso1,axs.flat):
    for aux in range(len(hs)):
        h = np.abs(ts[aux][0]-ts[aux][1])
        N = np.arange(0,len(ts[aux]))
        mat_parte1[num][aux] = metodo(lamda,h,N,x0)
        ax.plot(ts[aux],metodo(lamda,h,N,x0))
        #ax.set_ylim([-10,10])
        #Faltan los labels y titulos pero ya
        ax.set_title(f'{metodos[num]}')
    num  +=1
#Gráficacion
leg = [f'h = {pp}' for pp in hs]
axs[2].legend(leg,loc='center left', bbox_to_anchor=(1, 0.5))

<matplotlib.legend.Legend at 0x18771927b50>

<Figure size 640x480 with 0 Axes>

En este caso, el análisis de estabilidad nos indiga que h's van a converger a la función. Aquí se puede ver cómo principlamente en Euler forward la estabilidad dada por cada h influye mucho en que tan cercana se vuelve mi aproximación a la solución real

ks = [1,2,5,10]
hs = np.arange(0+0.1,1,0.1)
x0 = 10
def EF_caso2(h,k,xn):
    return xn + h*(-k*xn**2)
def EB_caso2(h,k,xn):
    a = h*k
    b = 1
    c= -xn
    return (-b+np.sqrt(b**2-4*a*c))/ (2*a)
def CN_caso2(h,k,xn):
    a = (h/2)*k 
    b = 1
    c= -k*xn+k*(h/2)*xn
    return (-b+np.sqrt(b**2-4*a*c)) / (2*a)

ts = [np.arange(0,10,i) for i in hs]


mat_parte1 = np.ones([3,len(hs)]) # 3 metodos* n hs * k * n tiempo
mat_parte1 = [list(mat_parte1[i]) for i in range(3)]
metodos_num_caso1 = [EF_caso2,EB_caso2,CN_caso2]

plt.figure()
fig,axs = plt.subplots(4,3,layout = 'constrained',figsize = (10,10))
nombre_metodo = ['EF','EB','CN']

num2 = 0
for k in ks:
    num = 0
    for metodo in metodos_num_caso1: 
        for aux in range(len(hs)):
            h = np.abs(ts[aux][0]-ts[aux][1])
            l_aux = [0]*len(ts[aux])
            #N = np.arange(0,len(ts[aux]))
            for i in range(len(ts[aux])):
                if i == 0:
                    l_aux[0] = x0
                else:
                    l_aux[i] = metodo(h,k,l_aux[i-1])
            mat_parte1[num][aux] = l_aux
            # mat_parte1[num][aux] = metodo(h,1,mat_parte1[num][aux])   
            axs[num2][num].plot(ts[aux],l_aux,label = f'h = {hs[aux]}')
            #Faltan los labels y titulos pero ya
        axs[num2][num].set_title(f'Metodo = {nombre_metodo[num]} , k = {k}')
        axs[num2][num].set_ylim([-5,10])
        if num2 == 0 and num == 2:
            axs[num2][num].legend(loc='center left', bbox_to_anchor=(1, 0.5))
        num  +=1
    num2 += 1

C:\Users\alanu\AppData\Local\Temp\ipykernel_27132\2543124031.py:5: RuntimeWarning: overflow encountered in double_scalars
  return xn + h*(-k*xn**2)

<Figure size 640x480 with 0 Axes>

F = [[0,1],[-1,0]]
x0 = [2,0]
hs = [0.1,0.5,1]
ts = [np.arange(0,10,i) for i in hs]
ns = [len(t) for t in ts]
l_x1 = [[[0]*n for n in ns] for _ in range(2)]
l_x2 = [[[0]*n for n in ns] for _ in range(2)]
def EF_caso3(F,x0,h,n):
    eigenvalues, eigenvectors = np.linalg.eig(F)
    E = eigenvectors
    D = np.diag(eigenvalues)
    E_inv = np.linalg.inv(eigenvectors)
    #A_diag = eigenvectors @ D @ np.linalg.inv(eigenvectors)
    I = np.identity(2)
    Xn =  np.matmul(  np.matmul(np.matmul(E,( np.linalg.matrix_power((I+h*D),n))),E_inv),x0)
    return Xn

def EB_caso3(F,x0,h,n):
    eigenvalues, eigenvectors = np.linalg.eig(F)
    E = eigenvectors
    D = np.diag(eigenvalues)
    E_inv = np.linalg.inv(eigenvectors)
    #A_diag = eigenvectors @ D @ np.linalg.inv(eigenvectors)
    I = np.identity(2)
    Xn = np.matmul( np.matmul(np.matmul(E,  np.linalg.matrix_power(np.linalg.inv(I+h*D), n)) ,E_inv)    ,x0)
    return Xn
#Solución análitica
def sol_analitica(t):
    return np.array([[2 * np.cos(t)], [2 * np.sin(t)]])
t_solr = np.arange(0,10,0.1)
aux_mat = sol_analitica(t_solr)
x1_analitico = aux_mat[0][0]
x2_analitico = aux_mat[1][0]


# Comparadas con la sol análitica
metodos = ['EF','EB']

plt.figure()
fig,axs = plt.subplots(3,2,layout = 'constrained',sharex= True)
num = 0
for metodo in [EF_caso3,EB_caso3]:
    num2 = 0
    for aux  in range(len(hs)):
        l_aux = [0]*ns[aux]
        #Paso
        for n in range(ns[aux]):
            l_aux[n] = metodo(F,x0,hs[aux],n)
        x1_num = np.real(np.array(l_aux).T[0])
        x2_num = np.real(np.array(l_aux).T[1])  
        l_x1[num][num2] = x1_num
        l_x2[num][num2] = x2_num

        axs[num2][0].plot(ts[aux],x1_num,label = f'h = {hs[aux]} ,m = {metodos[num]}')  
        axs[num2][1].plot(ts[aux],x2_num)
        axs[num2][0].plot(t_solr,x1_analitico,linewidth = 1.5,linestyle = 'dotted')  
        axs[num2][1].plot(t_solr,x2_analitico,linewidth = 1.5,linestyle = 'dotted')
        axs[num2][0].set_title(f'X1  h = {hs[aux]}')
        axs[num2][1].set_title(f'X2  h = {hs[aux]}')    
        num2 +=1

    num+=1
axs[0][1].legend(['Forward','Análitica','Backward'])
        

<matplotlib.legend.Legend at 0x1877244a890>

<Figure size 640x480 with 0 Axes>

plt.figure()
fig,axs = plt.subplots(1,2,layout ='constrained',figsize = (8,5))
num = 0
for ax in axs.flat:
    for j in range(3):
        ax.plot(l_x1[num][j],l_x2[num][j])
    num +=1
axs[0].plot(x1_analitico,x2_analitico,linewidth = 1.5,linestyle = 'dotted')
axs[1].plot(x1_analitico,x2_analitico,linewidth = 1.5,linestyle = 'dotted')
axs[0].set_title('Forward')
axs[1].set_title('Backward')
axs[1].legend(['h = 0.1','h = 0.5','h = 1','Análitico'])
axs[0].set_xlim([-5,5])
axs[0].set_ylim([-5,5])
plt.text(-0.4, 1.1, 'Espacio de Fase', horizontalalignment='center',
     verticalalignment='center', transform=ax.transAxes,fontsize = 13)

Text(-0.4, 1.1, 'Espacio de Fase')

<Figure size 640x480 with 0 Axes>

import numpy as np
import matplotlib.pyplot as plt

# Parámetros y configuración inicial
h = 0.1
r = 1
Ng = 15 
tiempo = np.arange(0, 10, h)
espacio = np.arange(0, Ng, r)
mat_parte3 = np.ones([len(tiempo), Ng])
def fun_inicial(x):
    s_asterisco = 5
    return (1 / np.cosh((x - s_asterisco) / 2)) ** 2

mat_parte3[0] = fun_inicial(espacio)

# Cálculo de la matriz
for j in range(len(tiempo) - 1):
    for i in range(Ng):
        i_p3 = (i + 3) % Ng
        i_p1 = (i + 1) % Ng
        i_m1 = (i - 1) % Ng
        i_m3 = (i - 3) % Ng
        if j == 0 and i == 0:
            print(i_p3)
            print(i_p1)
            print(i_m1)
            print(i_m3)
        mat_parte3[j + 1][i] = mat_parte3[j][i] + (h / (2 * r)) * (mat_parte3[j][i_p3] - mat_parte3[j][i_m1]) + (h / (8 * r**3)) * (mat_parte3[j][i_p3] - 3 * mat_parte3[j][i_p1] + 3 * mat_parte3[j][i_m1] - mat_parte3[j][i_m3])

# Graficar la matriz
fig, ax = plt.subplots()
cax = ax.imshow(mat_parte3, extent=[espacio[0], espacio[-1],tiempo[0], tiempo[-1]], cmap='gray',vmin = 0,vmax = 3) #, origin='lower'
cbar = fig.colorbar(cax, ax=ax)
ax.set_xlabel('Espacio')
ax.set_ylabel('Tiempo')
ax.set_title('Solitron')

plt.show()

3
1
14
12

from matplotlib.animation import FuncAnimation

x = range(Ng)
y = mat_parte3

# Inicializar la figura y los ejes
fig, ax = plt.subplots()
ax.set_xlim(0, 15)
ax.set_ylim(0, 1)

# Inicializar la línea que se actualizará en la animación
line, = ax.plot([], [], lw=2)

# Función de inicialización
def init():
    line.set_data([], [])
    return line,

# Función de actualización de la animación
def update(frame):
    x_data = range(Ng)
    y_data = mat_parte3[frame]
    line.set_data(x_data, y_data)
    return line,

# Crear la animación
ani = FuncAnimation(fig, update, frames=len(tiempo), init_func=init, blit=True)
ani.save( 's.gif', writer='imagemagick')
plt.show()

MovieWriter imagemagick unavailable; using Pillow instead.